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Abstract 

We present a theory of hnear optical constants based on a single-particle density matrix and im- 
plemented in an extension of the real-space muhiple scattering code FEFF. This approach avoids 
the need to compute wave-functions explicitly, and yields efhcient calculations for frequencies rang- 
ing from the IR to hard x-rays, and applicable to arbitrary aperiodic systems. Our approach is 
illustrated with calculations of optical properties and applications for several materials. 



I. INTRODUCTION 



This work is primarily concerned with theoretical calculations of optical constants, which 
are obtained from the long-wavelength limit g ^ of the dielectric function e{q,uj). These 
include the complex dielectric constant e{uj), the complex index of refraction, the energy- 
loss function, the photoabsorption coefficient, the photon scattering amplitude per atom, 
and the optical reflectivity. The ab initio calculation of these optical properties for arbi- 
trary materials has been a long-standing problem in condensed-matter physics.-i2i>^'^>^ Thus 
in practice, these properties are often estimated from atomic calculations or taken from 
tabulated sources .-i^i^i^ However, such tabulations are available only for a small number of 
materials over limited spectral ranges. Thus we aim to develop an efficient and widely appli- 
cable method covering a broad range of frequencies, thereby providing a practical alternative 
or complement to tabulated data. 

The theory of dielectric response of crystalline systems has been developed extensively 
over the past several decades,- following pioneering works of Nozieres and Pines,- Ehrenreich 
and Cohen,- Adler,- and Wiser.- These works developed the self-consistent field approach 
for the dielectric function within the time-dependent Hartree approximation, also known 
as the random phase approximation (RPA). Subsequently the theory has been extended to 
include exchange effects within the time-dependent density functional theory (TDDFT).— 
While ground-state DFT calculations are now routine, theoretical methods for accurate 
calculations of optical spectra are still not widely available. More elaborate theories have 
been developed that take into account quasi-particle effects and particle-hole interactions 
based on the Bethe-Salpeter equation (BSE),— li^iii^ but these are even more computationally 
demanding. 

In order to remedy this situation we have developed an efficient, real-space approach 
within the adiabatic local-density approximation that can be applied to arbitrary condensed 
systems over a broad range of frequencies from the visible to hard x-rays. Our approach 
is based on a density matrix formulation within an effective single-particle theory. This 
approach is a generalization of the real-space Green's function method implemented in the 
FEFF codes that includes both core- and valence-level spectra. Our work is intended to 
extend the capabilities and ease-of-use of FEFF to enable full spectrum output for general 
aperiodic systems with a quality roughly comparable to that in currently available tabulated 
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data.-^^^ This effort was begun by one of us using an atomic approximation for initial 
states.— This approximation is often adequate at high frequencies, but is unsatisfactory for 
optical and UV spectra. 

The remainder of this paper is arranged as follows. Sec. 11. describes the theoretical 
formalism behind our approach; Sec. 111. presents typical results for various optical constants 
for a number of materials; Sec. IV. discusses some additional applications and diagnostics, 
and Sec. V. presents a brief summary and conclusions. 



II. THEORY 



A. Density matrix theory of dielectric response 

We consider the macroscopic linear response of extended systems to an external electro- 
magnetic field of polarization e and frequency u 

Ve.t{t) = V,^t{u;)e^~'"'+''>' + cc, (1) 

where (5 is a positive infinitesimal corresponding to adiabatic turn-on of the perturbing po- 
tential. Throughout this work we use Hartree atomic units (/i = m = = = 1) unless 
otherwise specified. This perturbation polarizes the material, inducing a steady-state change 
Sn{'r,uj)e~'^'^'^ + cc in the microscopic electron density, which leads to a macroscopic polariza- 
tion P{uj)e~'^^^ + cc, representing the average screening dipole response of the electrons to 
the applied field. For simplicity of discussion, we assume that P has no component perpen- 
dicular to the applied electric field. This is the case for systems of cubic or higher symmetry 
in the g — > limit, but relaxing this restriction poses no computational difficulty. In this 
case one can define a scaler electric susceptibility x, and the dielectric function is — 

e(cj) = 1 + 47rx(cu) 

(2) 

where E is the electric field. Our calculations here make use of an effective single-particle 
microscopic theory in which the iV-electron state of the system at time t is described by a 
Slater determinant of time-dependent single-particle orbitals (j)i{t). Thus the state can be 
characterized by the single-particle density matrix p which is simply the projector onto the 
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orbit als: 

TV 

Pit) = Y,\m){m\- (3) 

1=1 

Their time evolution is governed by the time-dependent Schrodinger equation 

^j^\m) = Hm)) (4) 

for the time-dependent Kohn-Sham Hamihonian 

H = -iv^ + Kuc + Vh + Kc + + Kxt(t). (5) 

The terms in Eq. ((51) are respectively the kinetic energy, the electrostatic attraction to the 
nuclei Kiuc, the Hartree potential Vh, the ground-state exchange-correlation potential V^c, 
the dynamical contribution to the quasi-particle self-energy correction in the GW plasmon- 
pole approximation S^, and the time- dependent external potential Vext{t) of Equation (Q. 
Here and below, we suppress the position dependence of quantities when no confusion will 
result. The time evolution in Eq. (jl]) implies the Liouville equation^ for the density matrix 

»| = /fp-P^t. (6) 

In order to obtain the optical constants, we first linearize this equation with respect to the 
ground-state by decomposing the Hamiltonian and density matrix into their values in the 
ground-state and parts induced by Kxt 

H = Ho + H,{t) = Ho + V,xt{t) + l^nd(t) 

(7) 

P = po + pi- 

Hi consists of the external field and a term V^nd due to the response of the electrons. Second 
order terms, i.e., the products piHi and Hipi are discarded. We assume that the induced 
potential \4nd and hence Hi have the same time dependence as \4xt- With these assumptions, 
the time derivative in Eq. ([H]) becomes trivial and we can solve Eq. (jH]) for the induced density 
matrix in terms of the Kohn-Sham (KS) orbitals \(f)^) and eigenvalues Ei of the ground-state 
system, 

^-^ uj — [Lj — Li) + 10 

where fi = f{Ei) ^ 9{p — Ei) is the Fermi occupation number of state and p, is the 
Fermi level. The KS orbitals obey the unperturbed Schrodinger equation 

^1 mt)) = Ho m)) ■ (9) 



The induced electron density Sn{r, u) due to the perturbation Kxt is then given by 

5n{r,uj) = {f\pi(uj)\r) . (10) 

At this point it is convenient to introduce the bare and full susceptibilities whose local 
behavior is given by 

Sn{f,Lj) = {r\x°{ij)Hi\r) = (r|x(u;)Kxt|^ • (H) 

Typically, the bare response to an external perturbation is first computed from a single- 
particle (i.e. non-interacting) description of the ground state. The full response x of the 
system can be related to response x° of some non-interacting reference system. This proce- 
dure gives rise to the Dyson equation for x with an interaction kernel K 

X = X' + x'Kx = X°(l - ^X°)~'- (12) 

Methods for computing optical response that start from a single-particle description of the 
ground state can be classified by their approximations to the particle-hole interaction kernel 
K. The accuracy of the calculated macroscopic properties reflect that of the non-interacting 
response and the interaction kernel. Note, in particular, that one needs to find the frequency- 
dependent response of the non-interacting system, which involves different considerations 
than those for static, ground state properties (e.g., the ground state energy and density). 

In the crudest approximation K = 0: the resulting polarizability is that of the non- 
interacting reference system and local fields are neglected. In this case there is no screening, 
and the single-particle potential is the sum of the ground-state potential and An ob- 

vious deficiency of the non-interacting response is that the Coulomb field of the induced 
density is neglected. To address this deficiency Adler^ and Wiser-^ developed formally equiv- 
alent theories of the macroscopic dielectric response of periodic solids based on the RPA 
in which K is taken to be the bare Coulomb interaction. These theories were originally 
built on band structure calculations for periodic materials in the Hartree approximation, 
and Hartree local fields were included through the now termed Adler- Wiser formula. In this 
approach the operator inversion of Eq. (1121) is reformulated using the inverse of the micro- 
scopic dielectric matrix ecci^y k), which is then spatially averaged to give the macroscopic 
response e(a;) = limfe^o l/[e^^(^5 ^^)]o,o- However, the Adler- Wiser dielectric function is that 
of the Hartree system and has the deficiency that the underlying electronic wave function is 
not anti-symmetric under particle interchange. 
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Going beyond RPA thus requires additional exchange-correlation effects in K. There have 
been efforts along these lines of two types: those based on time-dependent formulations of 
density-functional theory (TDDFT), and those based on many body perturbation theory and 
the BSE. These approaches have been critically compared by Onida et al.^ By considering 
excited states from a quasi-particle viewpoint the interaction kernel can be decomposed 
into a direct term which is the Coulomb interaction between the quasi-particles and an 
exchange interaction , 

K = + K^. (13) 

Expanding Eq. (|T2l) in singly-excited (one electron, one hole) states and taking to be the 
Coulomb interaction screened by an effective (microscopic) dielectric function, yields a set of 
approximations referred to as the Bethe-Salpeter Equation (BSE). Various screening mod- 
els are used ranging from parametrized models (e.g. the Levine-Louie dielectric function) 
to independent-particle approximations such as the static RPA. BSE schemes can become 
computationally demanding since the inverse in Eq. (1121) must be dealt with in a product ba- 
sis which can be large. The differences between the independent-particle excitation energies 
and optical spectra and their interacting counterparts are referred to as excitonic effects. 
However, the non-locality of the exchange-correlation terms can be avoided by including 
exchange-correlation effects in K in terms of a density-functional /xc- Then the approach 
reduces to the TDDFT — where 

K{u) = v + Uc{^), /,,(^) = ^. (14) 

Consequently a local approximation to fxc leads to a local kernel (i.e., K depends only 
on the diagonal elements of the real-space single-particle density matrix). This locality 
imphes that Eq. ( fT2l) can be expanded in a single-particle basis, thus circumventing the 
need for particle-hole states needed for the BSE. The cost of this simplification is that direct 
information about the particle-hole interaction (e.g. exciton wave-functions) is only implicit. 
This makes it difficult to systematically improve on the local density approximation (LDA)i. 
Nevertheless, calculations in such TDLDA frameworks have been carried out for a variety of 
systems.— While the TDDFT has achieved good agreement with experiment for optical 
spectra in many cases, quantitative agreement at higher frequencies has been more elusive. 
Calculations with the BSE tend to be even more computationally limited. In addition these 
methods are built on various ground-state KS calculations, depending on the system. Each 
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approach can work well for a specific class of materials, but can lose accuracy or applicability 
for others. Also, the ground-state methods used were originally developed to calculate static 
properties and calculations of frequency-dependent (non-interacting) response can become 
cumbersome due to the need for large basis sets and special exchange-correlation functionals 
to describe unoccupied and excited states. 

The above difficulties have led us to consider a different approach with the goal of develop- 
ing a general method for calculations of optical response that can handle a variety of systems 
and spectral ranges. Our approach is based on an extension of real-space multiple scattering 
theory (RSMS) in terms of the one-particle density matrix. The RSMS approach is well 
suited to treat arbitrary aperiodic condensed-matter systems over a very broad frequency 
range (from the visible to hard x-rays). Indeed, this scattering-theoretic approach provides 
a superior basis for very high energy spectra where scattering is weak and the approach 
converges rapidly. Further the approach goes beyond the Born-Oppenheimer approximation 
and can include nuclear motion effects in terms of correlated Debye- Waller factors.— 

In this work, we present calculations within this RSMS approach using an independent 
quasi-particle approximation for the single particle states. Comparing Eq. ( ITOl) and ( |TT|1 
gives an expression for the bare response function or susceptibility 



Formally the imaginary part of the dielectric function is related to the full susceptibility 



where V is the volume of the system, and d = a ■ e^e is the transition operator between 
the incident photon of wave vector k and polarization e^. In practice the transition operator 
is replaced by the truncation to rank-one of its expansion into tensors developed by Grant,— 
which is equivalent to the dipole approximation. 

To evaluate Eq. (ITB]) for both optical and x-ray frequencies, we must first compute the 
response function x°(r, r cj). Formally Eq. (fT5|) can be expressed in terms of the single- 
particle Green's function as 




(15) 



byiS 




(16) 




(17) 



+ p{f',f,E)G'{f',f,E -u)dE. 
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Using the symmetries p{f, r E) = p{f', f, E) and G~{r, r E) = [G~^{f', f, E)]* on the real 
ii^-axis we can express the results entirely in terms of the one-particle density matrices p{E) 

Im y° /"-^^ 

—= p{f',r,E)p{r',f,E + Lj)dE. (18) 

^ Jep-w 

In this work we calculate these density matrices for energies ranging from the lowest occupied 
states to very high energies of order 100 KeV.— 

B. Multiple scattering Green's function 

Our calculations use an independent electron model in which each electron moves in an 
effective quasi-particle scattering potential V{f) which implicitly includes a dynamic self- 
energy correction T,d{E) to the ground state exchange and correlation potential. In this work 
T,fi{E) is calculated using the local GW plasmon-pole model of Hedin and Lundqvist.^^ The 
potential V{f) = J2n'^n{rn) + Vq is taken to be the self-consistent muffin-tin potential for a 
cluster of atoms at fixed locations Rn- Here r„ = r* — i?„ is the position relative to the n^^ 
atom, and Vq is a constant interstitial potential. Within RSMS theory, the Green's function 
for this potential can be written as a double angular momentum expansion 

G(f,r",E) = -2k[j2RLnirn)GLn,L'n'RL'n'ir'J 
LL' 

+ 6n,n' HLn{r>)RLn'{r<)] , (19) 

L 

where n and n' are the sites nearest r and r*' respectively, and r> (r<) is the larger (smaller) 
of the two position vectors. The terms in equation f|T9l) are the right-hand-side regular and 
irregular solutions Rin, H^n of the spherically symmetric single-site problems and their left- 
side counterparts -Rl„, -^Ln, the partial- wave phase shifts and the multiple scattering 
(MS) matrix GLn,L'n'- The wave functions are normalized so that in the interstitial region 
i?L„ coincides with li[/ij"'"e*'''" — h'^ e~^^'-"]/2i, and H^n coincides with Fi/i^e**^'". The bar 
for the left-sided solutions indicates that all factors except the Bessel functions are to be 
complex conjugated. Eq. (IT^ is rederived in the Appendix. As detailed in the Appendix, 
all these ingredients except the MS matrix can be found from the solution of a spherically 
symmetric single-particle quantum mechanics problem. The full MS matrix G for the system 
is found by numerical matrix inversion (e.g., with the LU or Lanczos algorithms in FEFF) 
with typical matrix dimensions of order 2 x 10^ or using the MS path expansion. 
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C. Relativistic basis 



To include relativistic effects such as spin-orbit coupling in our calculations properly it is 
necessary to recast the Green's function in terms of spinor solutions to the Dirac equation. 
In this context it is convenient to expand the spin-angular dependence of the one-electron 
states in the Pauli spinor-valued spin-orbit eigenfunctions which diagonalize both total and 
orbital angular momentum 



XKif-) = J2 Yp-''{f)rU\.m,-a,a\3,m\. (20) 



, 2 

' 2 



Here 0°^ is a Pauli spinor, K = {K,,mj) is a pair of relativistic angular momentum quantum 



numbers, and {I, s, mi, m^lj, rrij) is a Clebsh-Gordan coefficient. In this work as in Ref. 23| 



and [24], we have constructed the scattering matrix GLn,L'n' of Eq. ( 1191) using the scattering 
matrices ti calculated for the total angular momentum channel j = / + 1/2. This matrix is 
then transformed to the basis of spin-orbit eigenfunctions using Clebsch-Gordan coefficients. 
The central-site contribution Eq. (]A.17|) is constructed directly from numerical solutions of 



the central-site problem giving a total relativistic Green's function 

G{f,f',E) = -2k[j2HK n nn' 

KK' 

+ RKn{.^n)GKnK'n'RK'n'{.^'^) , (21) 

written in terms of right-hand (no bars) and left-hand (bars) solutions of the Dirac equation 
at energy E. These functions are 4-spinors which can be written in terms of the spin-orbit 
eigenfunctions: 

^ ... I ( P.K)x:'{r^) 

RKn{rn) = — 

^ ... if P.(r„)x™^^(r„) Y 

RKnirn) = — „ , , (22) 

where the T in Eq. ( l22l) denotes the transposed vector. The irregular solutions H, H take a 
similar form. These solutions are normalized by requiring the upper-component radial wave 
functions to coincide with [h'^e^^'^'^ — e^'^^''"]/2i (regular solution) or /i^e*"'"" (irregular 
solution). We are using the notation of Grant,— where the reader is referred for details 



regarding the numerical solutions P, Q appearing in Eq. (1221) . Tamura^ gives a relevant and 
illuminating discussion of solutions to the Dirac equation in spherical coordinates, although 
he treats a more general case using different notation. We have also transformed to a basis 
of real spherical harmonics to simplify calculations of the real-valued density matrices. 

D. Complex scattering potential 

The construction of the self-consistent muffin-tin scattering potential for the one-particle 
states is described elsewhere,— and we only briefly summarize the process here. First, a 
Dirac-Fock solver is used to calculate free-atomic potentials and densities which are then 
overlapped to obtain a starting point for the self-consistency loop. In this loop the one- 
particle Green's function for the full multiple scattering problem is calculated, from which 
a new electron density is calculated. Finally a new ground state muffin-tin potential is 
constructed within the LDA. The loop is iterated to self-consistency which typically takes 
about 10-20 iterations. Self-energy corrections are subsequently added for unoccupied states 
within the GW plasmon-pole approximation. 

E. Core state response 

At low energies (below the bottom of the valence band), the density matrix becomes 
sparse in energy, taking non-zero values only at isolated eigenvalues. In this regime, it is 
more computationally efficient to use orbitals to describe the electronic structure. Thus, we 
separate the single particle density matrix into two energy regions: the core region in which 
the atomic approximation is valid and the solid-state region where solid-state corrections 
are important, 



The core-valence separation energy Ecv is chosen to be an energy away from all KS eigen- 
values that separates the two regimes and is set by default to —40 eV, which is typically 
about 30 eV below the Fermi level. Above this energy p{E) is derived from the single- 
particle Green's function as described below. Note that in general there are occupied and 
unoccupied states above Ecv, but there are no unoccupied states below Ecv Similarly, the 




E < E, 



E> E, 



'CV • 



(23) 
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dielectric function e2(c<j)can be separated into contributions e2°'^^(c(j) and e^^cj) arising from 
transitions with core and valence initial states respectively. 

The core states are represented by single-particle atomic-like orbitals (p^. Here the index 
u = {n, i) denotes both a site index n and atomic level index i for the particular bound state 
at that site (e.g. Is, 2s, 2pi/2, etc.). We replace p{E) in Eq. (ITSl) for E < with 

PcoreiE)=J2pi\E), (24) 

V 

pi?(E)=0.(r>,(f')5(i?-e.). 
Thus we recover an expression equivalent to Fermi's golden rule for the absorption of light 

= — Im( i \d^G(uj + ei)d \ i )e(uj + ei-eF). 

V 

The initial core states \^y) and their associated eigenvalues are described accurately by 
Dirac-Fock atomic states for a single atomic configuration^"^-. For energies below i?cv ~ 
/i — 30 Ev < Vq the eigenf unctions of the central site problem are tightly bound to the 
central atom; their wave-functions decay rapidly as a function of the distance from the 
central site and can be taken to vanish in all cells except the central cell. This, along with 
the selection rules, limits the elements GxriK'n' (representing the final states) that contribute 
to absorption. For core initial states, the final state energy includes the inverse core-hole 
lifetime which broadens e'^. The calculation of the density matrix elements appearing 
in equation f|T6l) is handled differently depending on the photoelectron energy E = uj + e^,. 
For low-energy (less than ^ 50 eV+Vo) final states G is calculated by FMS just as in the 
calculation of e™^. At very high energies we again employ an atomic model and neglect 
scattering contributions (i.e. GkuK'ti' = in Eq. (12T]) ). At intermediate energies (50 eV 
+Vo < E < 1000 eV) we use efficient path filters^^ developed to treat EXAFS to find the 
dominant terms in the multiple scattering path expansion and sum these contributions to 
obtain the necessary GkuK'ti' elements, 

G = + G^TG^ + ■ ■ ■ . (26) 

The calculation of €2°'^'^ is accomplished by looping over the edges u with eigenvalues below 
i?cv For each edge we calculate €3 via FMS, path-expansion, and the atomic approximation 
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on appropriate energy grids. At this stage, correlated Debye- Waller factors can be included 
as in conventional XAS calculations using FEFF. 



F. Valence response 

Using the formal relation between the density matrix and the one-particle Green's func- 
tion p{E) = {—l/7c)lm.G{E) one obtains from Eq. (I2T1) 

p™^(f, r E)=J2 RKn{f^PKn,K'n'RK'nir'^'), (27) 
K,K' 

which is valid for r in cell n and f in cell n', where Rxni'f^) = XK{'f'n)RKn{i^n)- For real 
energies, the density matrix can be expressed entirely in terms of the regular solutions 
Rxn, and the irregular solutions do not enter. Below the Fermi level on the real energy 
axis, the density matrix is a rapidly varying function of energy. Away from the real axis, 
however, the behavior is much smoother. To both retain the separable form of Eq. (12 7p and 
the smoothness obtained by calculating the Green's function away from the real axis, we 
introduce a small broadening F and renormalize the regular solutions, so that the central 
atom density matrix gives the same density of states (DOS) in each Norman sphere as the 
actual broadened density matrix: 

i?,„(r; E) = A^n{E, F) Re [i?,„(r; E + tT)] , 

'"''Nrm 

(28) 



j,Nrm 








RKn{r\E) 


r'^dr = Im / 


/o 




Jo 



r'^dr 



X / r'\r'RUr<]E + iT)HUr>]E + iT). 
Jo 

This result is a key simplification in our approach. Here the Norman radius r^™ is defined as 
the radius of a neutral sphere centered on the rfi^ atom in the charge distribution formed by 
overlapping the charge distributions of the isolated atoms in their solid-state positions. The 
separable representation of the density matrix in Eq. fl271) permits a separation of the double 
spatial integral in Eq. (fT6|) into a product of two one-dimensional integrals. To complete the 
spatial integral in Eq. f|T6l) . we make the approximation that the spherical Norman cells n 
partition space and write the full integrals as sums of integrals over individual cells 

,(") 

JV 



dr — ^ X] / dr = ^ j rldrn j dVtn- (29) 
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The dipole matrix elements at each site n are defined as 

MIj,,{E,E')= [ dfRKn{r;E)dRKn{r;E')). (30) 

In the dipole approximation the matrix elements vanish except for transitions with j ' = 
j ± 1. Left (right) circularly polarized light only induces transitions with rrij' = m.j + 1 
{mj ' = rrij — 1). Thus the transition matrix M is sparse. Relaxing the dipole approximation 
is straightforward. Doing so introduces additional non-zero elements to M. With these 
conventions, the contribution to the spectrum from the response of the valence states (i.e. 
those occupied single-particle states with eigenvalues above Ecv) is given entirely in terms 
of density matrices and matrix elements, 

A l-Ep 

er'(^) = TT / dEY, Tr p„„,(E)M„,(E, E + uj) 

xpn'n{E + ij)Ml{E + u,E), (31) 

where and M„ are matrices in a truncated relativistic angular momentum K = [k, fi)- 
space. By symmetry, the sum over sites n in Eq. (!3T1) can be reduced to a sum over inequiv- 
alent sites in the solid. To compute e™^ we first solve the Dirac equation at each inequivalent 
site which yields T. Then Gmun' is found by inverting the full multiple scattering matrix, 
and matrix elements M are evaluated using the wave functions from the calculation of T. 
Finally, Eq. fl3Tl) is evaluated using trapezoid rule integration for the energy integrals. 

G. Spectrum construction 

With the response of both the valence band and the more tightly bound electrons calcu- 
lated, the contribution from each core edge is then interpolated onto a final output grid and 
combined with the other core edges and with the valence contribution: 



Y^et\uj)^et{uj\ (32) 



III. THEORETICAL OPTICAL CONSTANTS 

The examples presented here are primarily monatomic crystals (metals and insulators) 
with a single inequivalent site. However, the generalization to heterogeneous materials is 
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straightforward, and an example is also presented for AI2O3. Non-periodic materials can be 
treated by including enough sites to converge the spectrum. The calculations presented in 
this section used FMS matrices truncated at / = 3 and 147 atoms for all materials except 
diamond. The diamond calculation used / = 2 and 450 atoms. We include Diamond because 
it is a difficult case for our real-space method even though typical k-space calculations of 
Diamond (such as the plane wave psuedopotential calculation shown in Fig. [5]) use a unit 
cell containing only two atoms and can be less computationally demanding. All spectra 
were obtained by summing the contributions from 70 atoms. The response for the valence 
bands is obtained by calculating p™^ on a regular energy grid of 200 points. Then the dipole 
matrix elements M{E, E') are calculated for all pairs {E, E') with E below the Fermi level 
and E' above it. Eq. fl?I]) is then evaluated by matrix multiplication and simple numerical 
integration. To compute e^^cj) to high frequencies, we employ an atomic model of the 
valence bands based on average band energies and occupations calculated from p™'. The 
core state response is first calculated on a set of five 100 point frequency grids for each 
core initial state k, in the embedded-atom approximation. The FMS and path-expansion 
calculations are then carried out in cluster sizes of around 175 atoms on frequency grids of 
approximately 120 points. The contribution to €2 for each core initial state and the valence 
bands are then interpolated onto a large (5 x 10^ points) frequency grid which spans the full 
spectrum (e.g. 10~^ through 10^ eV) and serves as the final output grid. This grid has a 
higher density of points at low frequencies and around each core edge. 

A. Dielectric function: Imaginary part 

The fundamental quantity needed in our calculations of optical response is the imaginary 
part of the dielectric function €2(00) given by Eq. (|32l) . All other optical constants can be 
obtained in terms of e2{uj) as described below. As illustrative examples our density matrix 
calculations of €2(00) for Cu and Au are plotted in Fig. [1] compared to experiment. To 
demonstrate the effects of structural disorder on the dielectric response, we compare the 
imaginary part of the dielectric function for Diamond and amorphous Carbon in Fig. [32l 
Amorphous carbon structures were obtained with a "melt-and-quench" algorithm^^ using 
first principles molecular dynamics as implemented in the VASP package.-^^ These results, 
as well as the results presented below and calculations for other materials, are currently 
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FIG. 1: Calculated and experimental e2 for Cv^'^ (top), Au^^ (middle), and diamond^ and 
amorphous C'^^ (bottom). In the bottom panel the diamond curves have been shifted vertically for 
clarity. 

available in both graphical and tabular form on the FEFF website.— 
B. Dielectric function: Real part 

Owing to the analyticity of the dielectric response, the real and imaginary parts of the 
dielectric function are related by the Kramers-Kronig relatiort^ 

e{u;) = l + ^V Trf^'^^^^^. (33) 

Here V indicates the principal value of the integral. Since the denominator of the integrand in 
Eq. (I33l) has a pole aXui' = uj care must be taken when evaluating the transform numerically. 
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To evaluate the integral appearing in Eq. (133|) over the interval (cuj, cuj+i) between the i*^ 
and {i + 1)**^ grid points we find a linear approximation €2(0;') = muj' + b, which allows us 
to rewrite the Kramers- Kronig integral as follows: 

V / du = miooi^, - 00,) + (34) 

b — muj ^ fuji>i + uj\ b + muj ^ fuji^i — uj 
m H m ' 



UJi + UJ J 2 \ UJi ~ U! 

This expression is used to produce ei on the same output grid used for the imaginary part. 
The results of this procedure for diamond, Cu and Al203are plotted in Fig. O Even though 
the numerical transform Eq. ( |35l) is stable and accurate and (along with the calculated €2) 
completely determines ei via Eq. ( l33l) . we find that the real part of the dielectric function 
is more sensitive to errors and approximations than the imaginary part. 



C. Energy-loss 

With both real and imaginary parts of e{uj) one can easily obtain the energy loss function 

-lme-\uj)= ''^"^{^ (35) 

This is illustrated for Cu, AI2O3, and Au in Fig. [31 The loss function is proportional to the 
long- wavelength limit of the dynamic structure factor S{q,uj), which can be measured by 
inelastic scattering of either electrons in electron energy loss spectroscopy (EELS) or photons 
in non- resonant inelastic x-ray scattering (NRIXS). Calculations of the latter performed in a 
framework similar to ours have recently been reported by Soininen et. al. who only address 
the response of core electrons, but at finite q.^^ In contrast to ei we find that the loss function 
is less sensitive to errors and approximations in the density matrix than 62. Onida, et. al^., 
in an illuminating discussion of the differences between absorption and EELS experiments, 
have given an explanation of this observation in terms of the long-range part of the coulomb 
interaction. 

D. Index of refraction 

The complex index of refraction is simply the square root of the complex dielectric func- 
tion 

n{uj) + ik{uj) = e{ujy^\ (36) 
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FIG. 2: Calculated ei for diamond (top) and Cu (middle), and Al203(bottom) compared to ex- 
periment.— i2E 

Typical results for the real part of the index of refraction are given in Fig. |H 
E. Absorption coefficient 

The photon absorption coefficient fi{uj) is defined as the (natural) logarithm of the ratio 
of the incident and transmitted intensities for a photon beam across a thin sample, divided 
by the thickness. Theoretically fj^{uj) can be expressed in terms of the imaginary part of the 
index of refraction k{lj) 

/i(u;) = 2-k{uj). (37) 
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1000 



FIG. 3: Calculated energy-loss function (Eg. [35]) for Cu (top), AI2O3 (middle), and Au (bottom) 
compared to experiment.— 



Thus, fi{uj) is directly measurable with optical absorption experiments. Such experiments 
are currently performed to high accuracy using synchrotron light sources. We compare our 
calclulated results with experiment for several materials and with a calculation based on 
electronic structure calculated with ABINIT. This calculation was acomplished using the 
AI2NBSE package developed by Lawler, et. al.— which employs a BSE solver developed 
at NIST to generate optical spectra. The calculation shown excludes both local fields and 
excitonic effects and was generated using a regular grid of 8'^ fc-points to sample the Brillouin 
zone, 50 bands, and an energy cutoff of 30 Hartree for the plane wave basis. The C 1^ 
electrons were treated with a Troullier-Martins psuedopotential. For a sensible comparision, 
no gap corrections were included in either calculation. 
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FIG. 4: Calculated real index of refraction for Cu (top), AI2O3 (middle), and diamond (bottom) 
compared to experiment.— 

F. Reflectivity 

An important optical experiment for materials that can be prepared by vapor deposition 
methods is the measurement of the the reflectivity R defined as the ratio of the power re- 
flected from a planar face of a sample to the incident power. This quantity can be related to 
the dielectric response of the material by considering the boundary conditions satisfied by 
Maxwell's equations at the interface between the sample and vacuum. This procedure pro- 
duces the familiar Fresnel equations^^ relating the amplitudes of the transmitted (refracted) 
and reflected waves to the amplitude of the incident wave. As discussed by Stratton,— R 
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FIG. 5: Calculated absorption coefficient /x in inverse cm for Cu (top), Au (middle), and dia- 
mond (bottom) compared to experiment.—!^ The calculated diamond result is also compared to a 
reciprocal-space calculation. 



can be found by squaring the Fresnel equations. For example, for normal incidence 



[n{Lj) - if + n\u) 



(38) 



[n{Lu) + lY + k^{lu) 

The general expression for a lossy material (e2 7^ 0) and arbitrary angle of incidence is 
complex. However it is interesting to note that off normal incidence R{uj) has polarization 
dependence even for isotropic media. 



20 



12 
10 
8 
6 
4 
2 

-2 
-4 
-6 





! \ Theory 




;' \ Experiment 











10 



15 



20 



25 



30 




100 1000 10000 

(i)(eV) 



FIG. 6: Calculated real part of the anomalous atomic scattering factor for diamond (top), Cu 
(middle), and Au (bottom) compared to experiment.— '^^'^^ 

G. Photon scattering amplitude 



The Rayleigh forward scattering amplitude /{u) for photons can also be computed from 
the dielectric function^ 



UJ 



V 



e (co) - 1] . 



(39) 



Thus it is straightforward to calculate the x-ray scattering factors including anomalous terms 
using our RSMS approach f{q, u) = g{q, cj) + /*'*(g, u) + fi{u!)+if2{uj) in terms of /. Typical 
calculations of the real and imaginary parts of f{uj) are illustrated in Fig. 's [6] and [71 



21 



Theory 




Experiment 





15 
10 
5 



Theory 






Experiment 








>/ij LI / \\ i 
f'\\\\ j Y 















15 
10 
5 



1 — 1 1 — 1 ■ — ■— " 


-mq 1 

Theory 




Experiment 















1 10 100 1000 10000 

a)(eV) 

FIG. 7: Calculated imaginary part of the atomic scattering factor for Cu (top), Au (middle), and 
Al203(bottom) compared to experiment.— 

IV. APPLICATIONS AND DIAGNOSITCS 

A. Hamaker constant 

The Hamaker constant is the (real) function e{iuj) of a real frequency uj. For separation 
distances beyond the tunneling regime, the interaction between the tip and sample in an 
atomic force microscopy experiment is dominated by the van der Waals force, which can 
be calculated given the tip-sample geometry and the Hamaker constants of the tip and 
sample.— Using the analyticity of e in the upper half-plane, one can derive the following 
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FIG. 8: ncfT(w) calculated from the e2 sum rule for Mg using Eq. (j41[) . which assumes an asymptotic 
value of 12, the atomic number of Mg, in the w ^ oo limit. 



Kramers-Kronig type transform for the Hamaker constant 

/o 



e[iuj) = - duj — . 



(40) 



We evaluate Eq. fHU]) numerically in the same way we evaluate the Kramers-Kronig transform 
from €2 to ei, although away from uj = the integrand is regular. 

B. Sum rules 



Included in the output of our code are a few quantities useful for understanding the 
relationship between the underlying electronic structure and the frequency dependance of 
the optical constants. The /-sum rules for the imaginary parts of the dielectric function and 
the inverse dielectric function provide an important quantitative check of the calculation. 
We define the effective number of electrons per atom participating in transitions at frequency 



UJ 



V 



This quantity has the limii)^ 



lim Uesi^^ 



(41) 



(42) 



where Z is the number of electrons in the subsystem whose number density is N/V. The 
theory and calculations presented here are valid over a frequency range large enough to 
quantitatively evaluate the limit (jUj); missing or extra oscillator strength implies invalid 
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FIG. 9: JDOS/w^ foi p d (dashed line) and p s (dots) transitions and the calculated £2 of 
this work (solid line) for diamond vs. photon frequency in eV. 

approximations or unconverged calculations. Another check can be given by the index of 
refraction sumrule. 

POO 

[n{uj) -l]duj = 0. (43) 



C. JDOS 

As stated above, the selection rules constrain the angular momentum of final and initial 
states that can contribute to the absorption of light to a few channels (e.g. p —>■ d, s ^ p, 
etc.). The joint density of states (JDOS) correspondig to a certain dipole allowed chanel 
(/ —>■ I') is defined in terms of the normal /-projected DOS pf. 

/ pi{E)pv{E + u) dE, (44) 

J Ef—ui 

where the /-projected DOS is given in terms of the density matrix by 

Pi{E) = J2f df]YL{r)\'p{f,r,E). (45) 

m 

Neglecting energy dependence of the dipole matrix elements in the calculation of 62 gives 
a spectrum which is a sum of terms proportional to the JDOS/cj^ for the dipole allowed 
channels. We show this quantity for transitions from initial states with p character compared 
to the calculated €2 for Diamond in Fig. [91 
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V. CONCLUSIONS 



We have developed an efficient method for semi-quantitative ab initio calculations of op- 
tical constants over a broad spectrum, from the optical to x-ray frequencies. Our method, 
based on the one-particle density matrix, has been implemented in an extension of the RSGF 
approach in the FEFF codes which can be applied to general, aperiodic materials. We have 
illustrated the method here for a number of materials for which optical data are also avail- 
able including metals, insulators and aperiodic solids. Overall our results for the optical 
constants are semi-quantitative in the optical-UV range, but become much more quantita- 
tive for x-ray energies. Also their imaginary parts tend to be more accurate compared to 
experiment. This degree of accuracy is already adequate for many purposes, and especially 
for models which are not particularly sensitive to the detailed fine-structure in the spectra 
such as the calculation of screened coulomb potentials and van der Waals interactions. Fur- 
thermore many improvements are possible: i) It is desirable to include local field corrections 
as described above; ii) the muffin-tin approximation should be replaced with more accurate 
full potentials in each cell; iii) the extension to arbitrary momentum transfer q is often de- 
sirable. As noted above, the calculations can be done for any momentum transfer with only 
a modest increase in computational effort within our density-matrix formulation. In fact, 
the response of core states has already been extended to finite q by Soininen, et. al;— iv) for 
crystalline systems, it may be desirable and sensible to calculate the MS matrix in k-space, 
i.e., with periodic boundary conditions; and v) the treatment of the particle- hole interaction 
K currently only takes intra-atomic screening into account. 
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APPENDIX: REAL-SPACE MULTIPLE SCATTERING GREEN'S FUNCTION 

In this Appendix we describe the real-space Green's functions used in this work. Formally 
the Greens functions operator is given by 

G+{E) = [E - H + i5]-\ (A.l) 

where 5 is a positive infinitesimal. Expanding in the scattering potentials and free 
propagators yields the multiple scattering (MS) expansion 

G = GO + GVG = G° + G°TGO + ■ ■ ■ 

= [l-G-OTj-iGO (A.2) 
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Here we have introduced the local t-matrix tn = Vn + VnG^t to sum implicitly over all 
scatterings at a given site n, where (r|t„|r') = tn{r,r',E) vanishes outside a given cell n 
where w(r„)=0. 



1. Free propagator 

In position space the free propagator G^{E) is given by the FT, 

Below we evaluate this expression in terms of site-angular momentum scattering states \L, R) 
which diagonalize ti 



G'{r,r',E) = / ^J^^ J , (A.3) 



(A.4) 



jLifn) ={L,R\f)=t-'ji{krR)Yl{fn), 

where k = ^y2{E-Vo). 

In terms of spherical Bessel functions the free propagator is given everywhere by 

G^{r,r',E) = -2kY,yL{f)9i{r,r')Yl{f') (A.5) 

L 

= -2kY,ht{r>)jL{r<), (A.6) 

L 

where gi{r,r') = h^{kr^)ji{kr^) and h~j^{r) = i^h'l{kr)YL{r). This result can be obtained, 
e.g., from the FT using the identity exp(iA; ■ r) = 4TTJ]LjL{f)Yl{k) and carrying out the 
radial integrals in the complex /c-plane. Alternatively the same result follows from the 
inhomogeneous radial differential equation, where the prefactor is obtained from the Wron- 
skian 2/r'^W{ji, h^) = —2k. Here, as in the treatment of Rehr and Albers,^^ we have 
used the phase and normalization conventions of Messiah, with ji = {h'l — h'^)/2i and 
i^hi{x) = e^^ ci{l / ix) / X , ci is a polynomial of degree / with q(0) = 1. Also, for convenience, 
we have included the phase factors and in /i^ and Jl respectively, which do not change 
G", but simplify the asymptotic behavior. 

The expansion of the free propagator for points at different sites has the form of a matrix 
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product 

G%f,f',E) = Y.3L{rR)G\,L'R''jL{rR') 
= {f\LR) {LR\G%E)\L'R') {L'R'\f'). 

L,L' 

This follows directly from Eq. (1A.6P and the translation formulae for the spherical Hankel 
functions^ 

ht'{r'R) = J23LirR)G%^,,^,. (A.8) 

L 

Note the implicit factors of i'"' and z' in jL(rR) and jui^R) in this representation. In some 
works, e.g. that of Faulkner and Stocks^, these phase factors are included in the definition 
the propagator matrix elements. The above expression can be checked, e.g., by comparing 
i^h'l{kr) = j'likrjiji^' G'^jtmQ. Eq. (1A.7P can be derived, e.g., by expanding the exponen- 
tial product e*'^ (^'~'"') = e«^-(»"--f?)e-«^-(»"'-'f?')g«fc-(^--R') spherical Bessel functions, and then 
carrying out the integration over k. This procedure yields for the dimensionless propagator 
matrix elements: 

G1r,l'R' ^ = ^7rJ2{YLYL"\Y,.) hUkR"). (A.9) 

L" 

which depend explicitly on kR" = k{R — R'). The FEFF code uses dimensionless matrix 
elements Glj^,{kR) which have a separable representation^ 

jikR 

kR 



GlAkR) ^ G%,L'R' = Er^ArA,L', (A. 10) 



A 



pikR 



^ 4vr-^Qc;y;(i?)FiKi?), (^^ ^ oo), (A.ll) 

— * 

where Tx^i{kR) are generalized spherical harmonics. This can be obtained, for example, by 
substituting the asymptotic form of i''hi and and the completeness relation J2l YL{k)YL{R) = 
6{k-R). 

2. Full propagator 

Let us now evaluate the behavior of the full propagator G{f, r E) for r and r ' in different 
cells n and n' respectively. For this case the MS series can be viewed as a sequence of 
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scattering events consisting of all scatterings at site n followed by all sequences of scatterings 
not scattering at site n first or site n' last, followed by all scatterings at site n', 

Gnu' = [1 + G%-,]Gnn'[^ + t„,G°], (A.12) 

where the notation Gnn' refers to the propagator starting and ending in cells n and n' 
respectively, while Gnn' refers to those terms in the MS expansion with first scatterings 
at sites other than n and last scatterings at sites other than n'. This can be evaluated 
by substituting the representation of Eq. (1A.7I) into Eq. ( 1A.12I) and then re-expressing the 
terms in the site-angular momentum basis. Then Gnn' can be expressed in terms of the 
dimensionless full multiple scattering matrix elements Gin^Un' where 

jlirn) Gln^L'n' j li^n') 
L,L' (A.13) 

GLn,L'n' = [l " <5°T] ^ G°\^^^^,^, 

where <5°„^/„/ = — Snn')- The complementary delta-function in (5° ensures that 

G only includes initial scatterings from sites other than n and and final scatterings from 
sites other than n'. Next the terms on the left and the right sides of Eq. flA.12p can be 



expressed in terms of scattering states -RL„(r„). To see this note that matrix elements of 
the dimensionless t-matrices can be expressed in terms of phase shifts as 

{iL\in\3L') =iln^L,L' 

tin = e"^'" sin 6in. 



Then using the representation of in terms of Bessel functions in Eq. (IA.6p . one obtains 

{f][l + G%]\LR) = RLn{fn)e''^'' 

(A.15) 

= AMrn) + hl{rn)iin]YL{fn), (r„ > rf ), 

where Riir) = Rin{r)YL{f) . Asymptotically Rin{r) = [/i^e**^'" — /i;~e~*'^'"]/2i sm{kr — 
Itt/2 + 6in)/kr. For r < r™*, the radial states can be obtained from the regular solution to 
the radial equation, matched to the above result. Similarly one obtains {LR\{1 + tnG^)\f) = 
-RL(^n) exp(z(5i„). Note that the radial functions Rin{r) in the scattering states are real for 
real, nonnegative k, but are otherwise the analytic continuation to complex k. Combining 
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all these results in Eq. (1A.12I) then yields 

G{f,f',E) = -2k 



Rbnir T.)GLn,L'n'RL'n' ('^n')! 

LL' 

Glu^L'W = e''"^GLn,L'n'e''^'"' . (A.16) 



It is straightforward to show that this expression is equivalent to that of Faulkner and 
Stocks^3_ 

For r and f' at the same site n, G = G^ + G'^tnG'^ + Gn,n, where G is given by Eq. (lA.lSp . 
This yields 



G{r,r',E) = -2A; [ i/i„(f>)i?i(f<) 

L 

Ln,L'nRL'n (^i) , (A. 17) 

L,L' 

where Hiif) is the outgoing scattering state at site R which matches to i^e^^^^h^ {krn) for 
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